Download NHD Data

Here, we download NHDPlusV2.1, extract the Catchment polygons, and clip to the set of catchments that intersect HUC4 0303 (Cape Fear)

hu04 <- sf::read_sf("https://geoconnex.us/ref/hu04/0303") %>% st_transform(4326)
# download_nhdplusv2(outdir="./data/nhd/",
#                       url = paste0("https://s3.amazonaws.com/edap-nhdplus/NHDPlusV21/",
#     "Data/NationalData/NHDPlusV21_NationalData_Seamless", "_Geodatabase_Lower48_07.7z"),
#   progress = TRUE
# )
# 
# 
# path <- "./data/nhd/NHDPlusV21_National_Seamless_Flattened_Lower48.gdb"
# cat<- sf::st_read(path, layer = "Catchment") %>% st_transform(4326)
# 
sf_use_s2(FALSE)
## Spherical geometry (s2) switched off
# cat <- cat[hu04,]
# 
# write_sf(cat,"./data/staged/catchments.gpkg")

cat <- read_sf("./data/staged/catchments.gpkg")

#mapview(cat,alpha.regions=0, color="green") + mapview(hu04)

Download Census data

Here, we download census blocks, block groups, tracts, and counties from the Census TIGER/LINE files.

#  counties <- counties(state = 37) %>% st_transform(4326)
#  blocks <- blocks(state = 37, year = 2020) %>% st_transform(4326)
#  tracts <- tracts(state = 37, year = 2020) %>% st_transform(4326)
#  block_groups <- block_groups(state = 37, year = 2020) %>% st_transform(4326)
#  tracts_0303 <- tracts[hu04,]
#  blocks_0303 <- blocks[hu04,]
#  counties_0303 <- counties[hu04,]
#  block_groups_0303 <- block_groups[hu04,]
# # 
#  write_sf(blocks_0303,"./data/staged/blocks.gpkg")
#  write_sf(counties_0303,"./data/staged/counties.gpkg")
#  write_sf(block_groups_0303,"./data/staged/block_groups.gpkg")
#  write_sf(tracts_0303,"./data/staged/tracts.gpkg")

blocks <- read_sf("./data/staged/blocks.gpkg")
block_groups <- read_sf("./data/staged/block_groups.gpkg")
tracts <- read_sf("./data/staged/tracts.gpkg")
counties <- read_sf("./data/staged/counties.gpkg")

mapview(counties, col.regions="red", alpha.regions=0.2) + 
#  mapview(tracts, col.regions="orange", alpha.regions=0.2) + 
  mapview(block_groups, col.regions="green", alpha.regions=0.2)# + 
 # mapview(blocks, col.regions="blue", alpha.regions=0.2)

Intersection

We attempt to intersect blocks and catchments. First we ensure we are using valid spherical geometries, and calculate the areas of both. In the HUC4 0303, \(127,069\) of these result from \(50,379\) Census Blocks and \(15,357\) NHDPlusV2 Catchments. The histogram below visualizes the distribution of Census 2020 Block/ NHDPlusV2 Catchment intersection polygons according to the proportion of the area of the Block that is in the intersection polygon. This shows why using an intersection method is important, and not using simplifications such as assigning Blocks to Catchments based on centroid-in-polygon type methods. While \(44%\) of blocks in this sample are completely or \(>99%\) within one NHD Catchment, \(~32%\) of blocks are between \(20%\) and \(80%\) overlapping with 1 or more NHD Catchments.

sf_use_s2(FALSE)
cat <- st_make_valid(cat)
blocks <- st_make_valid(blocks)
cat$area_catchment_m2 <- st_area(cat)
blocks$area_block_m2 <- st_area(blocks)
cross <- st_intersection(cat,blocks)
## although coordinates are longitude/latitude, st_intersection assumes that they are planar
## Warning: attribute variables are assumed to be spatially constant throughout all
## geometries
cross$area_cross_m2 <- st_area(cross)
cross$prop_block_in_catchment <- as.numeric(cross$area_cross_m2/cross$area_block_m2)
cross$prop_catchment_in_block<- as.numeric(cross$area_cross_m2/cross$area_catchment_m2)

hist(cross$prop_block_in_catchment, breaks=10)

Areal Weigting Scheme

  1. ↩︎